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We introduce a high-performance implementation of a loosely coherent statistic sensitive to signals 
spanning a finite-dimensional manifold in parameter space. Results from full scale simulations on 
Gaussian noise are discussed, as well as implications for future searches for continuous gravitational 
waves. We demonstrate an improvement of more than an order of magnitude in analysis speed over 
previously available algorithms. As searches for continuous gravitational waves are computationally 
limited, the large speedup results in gain in sensitivity. 



PACS numbers: 07.05.Kf, 04.80.Nn, 95.55.Ym 
I. INTRODUCTION 



Loosely coherent algorithms [T] detect families of noise- 
dominated signals. The immediate application is the 
search for continuous gravitational waves. The large pa- 
rameter space volume to be searched and the data vol- 
umes necessary to search for exceedingly weak signlas 
make our searches computationally limited. Hence im- 
provements in efficiency of our algorithms directly affect 
the sensitivity of our results. As we are still waiting for 
the first detection, we cannot rely on a natural source 
to verify correctness of the detector and search pipelines. 
Our algorithms must be robust to possible imperfections 
of the detector, to faults in understanding of gravitation 
or even to bugs in the search programs. It also helps to 
be sensitive to a wide family of signals, in case the loud- 
est source is not a perfect sine wave, which can result, 
for example, from a companion object. 

The first implementation of a loosely coherent search 
[T] was designed for signals with a large amount of phase 
deviation over 30-minute interval from a perfect sine 
wave. This provided the gain in sensitivity needed for 
follow up of outliers seen in the full dataset of the LIGO 
detector's fifth science run S5 while preserving the 
robustness of the underlying, semi-coherent PowerFlux 
algorithm [3H5]. 

In this paper we explore the other end of the spectrum 
- an algorithm sensitive to coherent signals described by 
a small number of parameters, such as frequency or sky 
position. A number of coherent codes have been devel- 
oped previously, in particular [BHT3] . What is different in 
our approach is that, unlike previous algorithms, our sky 
templates are "thick" , sweeping small patches of param- 
eter spaces. In particular, individual signals taken from 
the middle of nearby patches do not have a high degree 
of overlap. This property, together with careful attention 
to implementation particulars, provides for a very high 
performance coherent code. 



II. A SIMPLIFIED ALGORITHM 

Suppose that we are interested in a family of signals 
described by the formula 

s(t,a) =^e2'^*''(*+-(*'"» (1) 

where v is the signal frequency and S(t, a) has a smooth 
dependence on time t and multidimensional parameter 
set a. 

Our input data f{t) — s{t, a) +^{t) consists of one sig- 
nal from this family and, ideally, uncorrelated Gaussian 
noise of variance cr^(t): 

mm) = ^'mt-t') (2) 

If we knew the frequency v and parameter set a, we 
could form a matched filter that would return the ampli- 
tude of our signal: 

Here W is the total weight: 

J to W 

If the signal parameters are not known, one can con- 
struct a bank of waveforms s(t) and evaluate the integral 
separately for each template. This is, of course, compu- 
tationally expensive. 

One way to gain a large speedup is to introduce a new 
time variable t' that straightens out our signal into a sine 
wave: 

t' = t + E{t,a) (4) 

One then resamples [T^ the input data /(t) to be equally 
spaced in the new variable t' and uses a Fourier transform 
to compute amplitudes for a range of possible signal fre- 
quencies. It has also been proposed by B.F.Schutz that 
a change to nearby value a' can be accomplished with a 
kernel [16]. As far as we know this method of "stepping 
around in the sky" does not yet have an implementation. 
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A simpler approach that combines the best features 
of both resampling and stepping around in the sky is to 
consider the following function of the input data: 



F(A;i/o,ao) = 



1 



/(^) -27TWo(t+S{t,ao)) -2TTi\t 



dt 
(5) 

which is easily computed with a fast Fourier transform. 
For A = it returns an amplitude estimate of the signal 
with parameters (I'oi'^o)- The values of F for X ^ 
carry slightly distorted information on nearby templates, 
which can be used to compute estimates of their signal 
amplitude with a convolution: 



1 



dt 



2U 



/We 



-2TviiJo(t+S{t,ao)) 



to. 



^g— 27r'i{i/— i/o)^— 27ri(i/ — i/o)!E;(f:,ao) — 27r2^'(H(t,a) — H(i,ao))^^ 



tl 



dtd^ 
(6) 

The reader will notice that the last term is not quite 
the ordinary convolution - one of the convolved terms 
has a (slow) dependence on the convolution parameter. 
We call this a "pseudo-convolution" operator. We distin- 
guish this case from the more general notion of integral 
operator, because in practical computation we do not 
have to update the slowly changing convolution with ev- 
ery data sample and the computational requirements of 
the operator are equivalent to the computational require- 
ments of plain convolution. 

In our case the pseudo-convolutions are of the form 



A{v, a) 



F{v-iyo-^;vo,ao) 



-2TTi^t — U{t,y,a) 



dtdfi 

. 

where U{t, v, a) is a phase mismatch function describing 
difference in phase evolution between nearby templates. 
For smooth U{t, v, a) the convolution operator is close 
to (5-function. In practical computation, using discrete 
Fourier transform, this means that our convolution can 
be approximated with an FIR filter that has small num- 
ber of terms. 

It is crucial to control the number of terms in the con- 
volutions. There are two ways to achieve that, aside from 
simply using small increments in the parameters with a 
corresponding increase in the number of templates. 

First, we can subtract a linear term from the argument 
of the exponent so that it has the same value at both 
ends of the segment [to,ti]. The linear term is analogous 
to a Doppler shift correction and results in relabeling of 
frequency parameter v. 

Secondly, pseudo-convolutions have small or null com- 
mutators. One can then apply methods of linear alge- 
bra to change from initial set of operators (usually corre- 



sponding to individual parameters) to a set with progres- 
sively fewer convolution terms. The operators with the 
smallest number of terms are used in the innermost com- 
putational loop thus determining overall performance of 
the code. 

We implement these techniques by representing 
J7(i, V, a) as a sum of a linear term and Fourier series 
with coefficients linear in v. 



^ _ U{ti,v,a) -U{to,v,a) 
U{t,i^,a) ~ (t-to) + 



tl — to 

oo 

+ J2 {Maf + ul{a)i 



Jl-nikt 



(8) 



fc— — oo 



In most cases the equality is exact and only one of vP^ or 
u\ is non-zero. The set of pseudo-convolution operators 
can be transformed into a more convenient basis by min- 
imizing coefficients in the series, which is then converted 
to its Fourier transform with the help of the Jacobi- Anger 
identity: 



E 



z e 



(9) 



Only a small number of terms are usually needed and 
recomputation is done rarely. The simulations presented 
in section [TV] were carried out with innermost loops that 
used convolutions with only 11 terms - a number chosen 
to take advantage of vectorized arithmetic on modern 
CPUs. 



III. IMPLEMENTATION DETAILS 

While a fast engine to compute coherent sums is es- 
sential for our search code, it is only part of a whole. In 
particular, after computing coherent sums one needs to 
derive statistics such as maximum SNR or upper limit 
which can be expensive to compute. For example, if one 
were to use a rank-based method (which is nicely robust), 
it would require sorting the computed sum which has a 
scaling of N log N - same as a fast Fourier transform, and 
much slower than a convolution. 

It is usually not practical to analyze the entire band 
of interest in one go, but rather one splits it into fre- 
quency bands of 1 Hz or smaller. The amount of loaded 
data can be greatly reduced by precomputing short dis- 
crete Fourier transforms (SFTs) of duration commensu- 
rate with the region of interest. It is convenient to have 
the SFT length be short enough that the signal frequency 
can be assumed to be stationary. 



A. Polarization analysis 

Continuous gravitational waves have a more compli- 
cated form than is given by equation [l] - there are two 
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polarizations with detected strengths that vary with the 
orientation of the detector. 

The following analysis is similar to one found in 
[51 [T71 [THj; we prefer, however, to reduce the four real 
parameters to two complex parameters that have a sym- 
metric role. We also derive a convenient equation for 
surfaces of constant Hq. 

We start by assuming that our signal consists of two 
polarizations: 



h'_^ = cos{ujt + 
h'^ = Ax sm{ujt + I 



(10) 



A generic pulsar signal can be represented as = 
/lo (l + cos^(t)) /2, Ax = /locos(t), with ho = A^ + 



A\ - and cos(t) ^ A^j [Aj^ + ^^ J^^ - A\ 

We will assume that demodulation is performed for a 
fixed frame of plus and cross polarizations rotated at an 
angle /3. In this coordinate system we have: 

/i+ = cos(a;t + 0) cos(e) — A^ sin(wt + sin(e) 
hx = A+ cos{ujt + (j>) sin(e) + Ax sin{ujt + cj)) cos(e) 

(11) 

where we introduced e = 2(-i/' — /?). 

The signal amplitude in SFT bin corresponding to fre- 
quency uj is then 

{F+h++Fxhx)e~''^*dt = 

= \e"l' {F+{A+cos{e) + iAxsMe))+ (12) 

+Fx {A+ sin(e) - iAx cos(e))) 
= Fj^wi + FxW2 

where we introduced complex amplitude parameters 

wi = le*'*(A+ cos(e) + iAx sin(e)) 



W2 = 2^ (^+ sin(e) — iAx cos(e)) 



(13) 



The complex amplitude parameters Wi and W2 are alge- 
braically symmetric: 

b = \w^\' + \w2\'^liAl+Al)^ 

= — /i^(l -|-6cos2(t) +cos''(i)) 
1 1 

a = Im(wiU'2) = -A^Ax = -/iq(1 + cos^(t)) cos(t) 
4 8 

. 

and are otherwise unconstrained. They have a simple 
relation to more familiar real amplitude parameters A'^ 

m- 



wi = - iA^ 
W2 = A^-iA"^ 



(15) 



One easily finds the following equation of constant ha 



y/\wi - iw2\ 



ho 



(16) 



the solutions of which form a singular surface enclosing 
a non-convex solid. This complicated form is responsi- 
ble for differences between PowerFlux style upper limits, 
which are always limited by sensitivity to linearly polar- 
ized signals, and SNR statistics, the outliers of which can 
have arbitrary polarization. A related issue is the differ- 
ence between the J^-statistic and the i3-statistic |13) . 
It is easy to compute generators for rotations in (p and 

e: 



d_ 
84) 

d f wi 
de \ W2 



Wi 
W2 



Wi 
W2 



-1 

1 



Wl 
W-? 



(17) 



(18) 



This shows that the surface oi ho — 1 is obtained by 
revolving the parabola 



along (j) and e. 



W[ = 1(1+COS2(0) 
W2 = — |cOs(t) 



B. Coherent sum 



(19) 



In a coherent analysis we have data for many SFTs 
{zi}^-^ which contain our signal mixed with instrumental 
noise 

z, = iF+iU)wi + Fx {U)w2) e'*(*') + e< (20) 
and we construct a weighted sum 



N 



- 'F+{U)w[+FxiU)w'2 



(21) 



which estimates signal amplitude. Here <&(ti) describes 
some assumed phase evolution due to changes in the 
source or detector, ai are weights satisfying X^iLi Q^i — Ij 
and w[ and are computed for polarization and phase 
of our signal, but assuming /iq = 1, i.e. they satisfy 



W\ + IWn + \ \W\ — IW. 



(22) 



There are many ways to compute "optimal" weights 
ai , in particular [6l 113) . Here we use the variance of Z 
as the optimality measure. Assuming are independent 
Gaussian variables with zero mean, we compute: 



N 



Var (Z) = V a? 



Var (C.) 



^ '\F+iU)w[+FxiU)w'2\^ 
One easily finds that Var (Z) is minimized for 
1 \F+{ti)w[+Fx{U)w'2\^ 



^{W[,W'2) 



Var (6) 



(23) 



(24) 
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where 21 is the normahzation weight: 



N 



|F+(<,K + Fx(t.K|' 



Var(C,) 

Substituting the optimal coefficients a^, we compute 



(25) 



JV 



Z{w[,w'2] 



and 



Zie 



i^(u) F+iU)w[+F^{U)w'2 



Var 



Var(ZK,^)) 



(26) 



(27) 



We see that the total weight ^{w'l , W2) can be interpreted 
as a measure of the amount of data used to compute 
Z{'w'i,W2), as in the case of stationary input data it is 
proportional to the number of independent SFTs. 

We now note that 2l(wi , W2) depends quadratically on 
coefficients Wi and on the estimate of variance of the data 
and can thus be computed once for all templates with 
similar detector response and Fx ■ The unnormalized 
coherent sum 21^ is linear in Wi, and these coefficients 
are exactly the kind of sum that we learned to compute 
in section |TT1 

Our computation then results in the following repre- 
sentation of Z: 



Z{n) = 



(28) 



Here Yii, Y12 and Y22 are coefficients that determine the 
total weight of the coherent sum and Xi (n) and X2 [n) 
are two arrays (in frequency bin n) of coherent sums: 



F+{Uf 
tVar(e.; 



Yr. 



Y22 — 



N 

E 



F+{U)Fx{U) 



Var(C. 



JV 

E 



- Var (e. 



(29) 



Xi 



X2 = 



N 

E 

1=1 

N 

E^ 



Var(C.) 
Var (CO 



For convenience, we tabulate a few useful expressions 
using these coefficients: 



^{W[,W'2) 
Z{w[,w'2) 



= FiihiP + 2yi2Re(u.iu;2') + ^22^1 



SNR(u'i,^) = 
SNRfl,(wi,W2) = 



Xiw[ + X2W2 
2l(w'i, w^) 

^{w[,w'2) 

(Re iX,w[ + X2w'2)f 

^{W[,W'2) 



(30) 



C. Efficient computation of coherent sum statistics 

We now need to reduce the data to the SNR or 
some other statistic of the loudest outlier. The simplest 
method, employed by PowerFlux [3]-[5] and the large-S 
loosely coherent search [Tj is to scan different values of 
wi and W2 looking for a maximum. The elements of our 
coherent sums are computed, however, with ~ 22 com- 
plex multiplications for each frequency bin, and even a 
modest grid of polarization parameters dominates com- 
putation. 

An approach taken in is to analytically maximize 
SNR over wi and 'W2. A new *B-statistic was introduced 
in 13! that was shown to have a better physically moti- 
vated prior. 

There is an easy and elegant way to compute all of 
these statistics with minimal cost. 

First we note, that our statistics are monotonic in sig- 
nal strength; they just disagree as to which signal pa- 
rameters are given preference. Once the signal strength 
is fixed (in any statistic), though, the rest of the pa- 
rameters form a bounded manifold - and other statistics 
achieve a maximum and minimum value on it. We can 
thus infer an estimate of another statistic from knowing 
the maximum of some convenient, easy-to-compute mea- 
sure of signal strength. 

As a toy example, assume that our statistics vary by at 
most a factor of 4 for signals of the same power -|- 
IX2P and that our signal array consists of w 300000 
elements. The average maximum of 300000 of mean=l, 
exponentially distributed samples is 12. But 95% of these 
numbers are below 3 = 12/4. Thus, to find the maximum 
of a more complicated statistics we only need to examine 
5% of the samples, providing a factor of 20 speed up. If 
a signal is present our maximum is even higher excluding 
a larger amount of data. 

In a practical implementation it is better to use ad- 
justed power (see table |l]), which results in a less-than-4 
worst-case scale factor for upper limit statistics (evalu- 
ated for strong signals) and close to unity factors for plain 
SNR and 5^-statistics as can be seen on figure [T] The 
strong dependence on declination is due to the antenna 
pattern of the detector. The fraction of templates for 
which we computed upper limit statistic during a Gaus- 
sian noise simulation run (discussed in more detail in the 
following section) is shown in figure [2j 
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Statistic Formula 



Raw power jXij^ + |X2l^ 
Adjusted power Y22 \X^\' + Vii {Xif 

CMC \Xlw'i + X2W2f 

SNR max 



\X^W[+X2w',f \X^W[+X2W'^\ 



Upper limit inax ^/ - — ^^] ^ n-/ + 2 , ,..,v^/9 + 
J-stat 



2 



y22|Xip - 2yi2Re(XiX2) + Kii|X,'2 



^11^22 - 



TABLE I: Statistics functions. The variables w'l and w'2 are constrained by equation |22| 



ratio_SNR o ratio_UL_circ x 
ratio_UL + ratio_F_stat <> 




-1.5 -1.0 -0.5 0.0 0.5 1.0 1.5 -1.5 -1.0 -0.5 0.0 0.5 1.0 1.5 

Declination, rad Declination, rad 



FIG. 1: Maximum variance of statistics for signals of constant 
adjusted power. The numbers come from simulations using 
Gaussian noise and describe the ratio between the maximum 
and minimum statistic values for signal of constant adjusted 
power. Top curve - upper limit assuming circular polarization. 
Next curve below is upper limit statistic, followed by SNR and 
5-statistic (color online). 



IV. PERFORMANCE AND VALIDATION 

An initial implementation of the ideas discussed above 
has been completed. It has a number of simplifications 
compared to an eventual production program - the input 
data is assumed contiguous, one spindown is analyzed at 
a time, and there is no support for higher-order source 
frequency evolution parameters. 

We have performed Monte-Carlo runs of 1000 injec- 



FIG. 2: Fraction of templates resulting in statistic calls. The 
underlying data was pure Gaussian noise (color online). 



tions each into Gaussian data spanning four million sec- 
onds (approximately 1.5 months), assuming a detector 
located at LIGO Hanford observatory. The injection 
sky locations and source orientation were uniformly dis- 
tributed. The spindown parameter was set to be zero. In 
the first run, we injected signals of various strength (fig- 
ure [s]) to test signal detection and upper limit estimation. 
The second run had identical parameters and noise dis- 
tribution, but the signal strengths were set to 0. This 
provided upper limits on pure noise alone (for relative 
comparison) as well as worst-case computational perfor- 
mance, as the presence of strong signals makes compu- 
tation of statistics faster. The injection frequency was 
uniformly distributed in ±0.084 Hz interval centered on 
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400 Hz. 

Each separate injection run included a search over a 
1-arcminute disk around a nominal sky location that 
was obtained from rounding the true injection locations; 
hence the actual injection locations were uniformly dis- 
tributed in relation to the sampled grid. 




-24.6 -24.4 -24.2 

Injected strain 



FIG. 3: Upper limit versus injected strain (color online). 



The upper limit plot in Figure [3] shows that upper lim- 
its are consistently above the injected signal values. A 
gap between the reconstructed points and the red line 
marking injected values is due to a conservative correc- 
tion for Hann-windowed input data. 

When signals rise above background their frequencies 
are well localized (figure [4|. We use this localization as a 
criterion for detection: a signal is considered to be found 
if its frequency is within 1x10^^ Hz of true value. This 
corresponds to false alarm ratio of ~ 6xl0~^. 

Figure [5] compares the efficiencies of various detection 
statistics. The SNR and iJ-statistics are mathematically 
equivalent; the only difference is that the SNR is com- 
puted by iterating over a grid of parameters w[ and 
while the S'-statistic has a much more efficient closed 
form. The slower SNR algorithm was used as a bridge 
to an implementation of the S-statistic, and will also be 
useful for future implementation of S > loosely coherent 
statistics. 

As seen from the plot we observe no difference in per- 
formance between S^-statistic and QS-statistic. We believe 
this is due to two factors: 

• First, even for pure noise, maximizing over 2x10^ 
bins results in high SNR values where the difference 
between g'-statistic and *B-statistic is smaller. 



nnxssmxsco ogqogd 



0.0 0.5 1.0 1.5 2.0 2.; 

Injected strain/Upper limit 
FIG. 4: Frequency reconstruction (color online). 



SNR 

B-statlstIc 



F-statistIc X 



^ 60 




1.0 1.5 2.0 

Injected strain/Upper limit 



FIG. 5: Statistics efficiency versus strength of injected signal, 
relative to established upper limit (color online). 



• Secondly, the area searched is small which sig- 
nificantly reduces infiuence of the weight term 
^{w'i,w'2) that appears in the *B-statistic 

It might be possible to take advantage of the improved 
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performance of the S-statistic at low SNR by analyzing 
multi-detector data with consistency cuts to bring down 
maximum SNR. 

The computational performance of our code is shown 
in Figure [6j The y-axis is in units of cycles per tem- 
plate, with each template computed once per frequency, 
sky position and spindown while sampling all possible 
alignments of the source. All statistics (SNR, upper 
limit, circular upper limit, 5^-statistic and *B-statistic) 
were computed during the run. The underlying data was 
purely Gaussian. The simulations were run on a cluster 
of 2.3 GHz AMD processors. 

Our worst-case performance is below 1500 cycles which 
favorably compares with performance of resampling [141 
[20] that is estimated to be « 20000 cycles per template. 
The cpu utilization by different parts of the algorithms 
is shown in table |ll] Only a third of the cycles is at- 
tributable to computation of the convolution, while at 
least 46% is spent gathering statistics, leaving much room 
for further improvement. 



CPU fraction 


Code description 


36% 


Computation of upper limit 


16% 


Computation of 11-term convolution 


10% 


Computation of terms of convolutions 


8% 


General statistics function 


7% 


Computation of 63-term pseudo-convolution 


2% 


Computation of logarithm 


Balance of 21% 


Spread out over many parts of the algorithm 



TABLE II: CPU cycles spent in different parts of the algo- 
rithm for sky position with declination of 1.0 radian. 



There are several directions for further improvement: 

• Making use of coherent or loosely coherent com- 
bination of data from two interferometers should 
improve sensitivity and provide better rejection of 
detector artifacts. This might also be an area where 
the *B-statistic will show its strength. 




-1.5 -1.0 -0.5 0.0 0.5 1.0 

Declination, rad 



FIG. 6: CPU cycles spent per template as a function of dec- 
lination of the injected signal (color online). 



• It is necessary to derive more efficient alternatives 
for the computation of upper limit. 

• It would be desirable to extend the algorithm to 
search over frequency evolution parameters to be 
able to detect binary systems with weak modula- 
tion. 

• The fine granularity of input data can be used to 
avoid high-intensity glitches, by excluding contam- 
inated SFTs. 
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V. CONCLUSIONS 

We have described an implementation of the loosely 
coherent statistic that searches over families of ideal con- 
tinuous gravitational signals. The performance of the al- 
gorithm is more than an order of magnitude better than 
previously published algorithms, opening the way for ex- 
ploring wider parameter spaces. 



Appendix A: Efficient computation of ©-statistic 

*B-statistic was introduced in [T3] as a Bayesian alter- 
native to 5^-statistic. It was shown that the 5^-statistic is 
equivalent to a Bayesian statistic with a prior that favors 
linearly polarized signals and high signal strength. In 
contrast, the 58-statistic is isotropic in spin orientation. 
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The statistic starts with the likehhood function: 



easy identities: 



_ ^ReiwiXi+W2X2)-^{Yii\wi\^+2Yi2ReiwiW2)+Y22\w2\^) 

(Al) 

Instead of maximizing it, which is the approach of the 
iJ-statistic we compute the integral: 



B{x) = / dAC{x,A) = dw' dh- 

Jh<ho J Jo 

.ghRe{iD[X^+iB'2X2)-^h^{Yll\w[\^+2Yl2'Re(w[^B'2)+Y22\w'2\^) 

(A2) 

The measures dA and dw' are chosen to be uniform in 
parameters 0, tjj and cos(t). The integral with respect to 
h is not normalized, which allows to set its upper limit 
infinite. This effectively changes B{x) to be in units of 
strain. It would be interesting to explore the possibility 
of deriving an upper limit estimator based on the same 
principles as B{x). 

The integral with respect to h can be shown to be a 
function of SNRfl(wi, and total weight ^{w[,w'2): 

POO -1 2 pa/^/b 

/ dxe"''-'"''/^ = -=e^ / dxe-"'/2 ^^3) 
Jo vb J-00 



B{x) 



dw' , =^ / 6-=" /^dx 



(A4) 

This still leaves a three dimensional integral to carry out 
which is undesirable inside the inner loop. We note that 
2t(wi,u;2) does not depend on phase ^, while 

SNRrK, w^) = cos(,^)SNR(w;, w^) (A5) 

We can thus represent our statistic as 

B{x)= [ dw' \ e(^SNRK,u;^)) (A6) 

J0=O ^y^{w'l,w'2) \^ J 



where 6(a;) is defined as 



e(x) 



1 f'^^ gCosl 



2-K 



{4>)x /•cos{cf>)x 



2-K 



dsd(t> (A7) 



<d{x) = Q{-x) 
G(0) = i 



(A8) 



It is also easy to compute approximations for small and 
large x: 



e{x) = -^(i + o{i/x)) 



(A9) 



Armed with these relations, we can spend some time in 
numerical experimentation and arrive at the following 
approximation: 



^ , , exp (\/0.25 + a;2) + a2x'^ + aiX'^ + x^ 

\y\X) ^ 

(47r2a;2 + ISe^)^/^ &o + h^x'^ + 642;^ + a;^ 

(AlO) 

which has a 0.05% error over the entire range with the 
following values of constants: 



a2 
a4 
bo 
b2 
bi 



7.7199014890487 
19.0337266315871 

5.2017224760755 
7.7201854234519 
21.1533518190664 
4.2818853782852 



(All) 



Now we can compute !B-statistic with a simple sum over 

a uniform grid in ■0 and cos(/,), at a cost within an order 
of magnitude of computing the SNR statistic. 



We can now study Q{x) as a new special function and find 
a means to compute it efficiently. We have the following 
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